Dressed Counterions: Strong Electrostatic Coupling in the Presence of Salt 
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We reformulate the theory of strong electrostatic coupling in order to describe an asymmetric elec- 
trolyte solution of monovalent salt ions and polyvalent counterions using field-theoretical techniques 
and Monte-Carlo simulations. The theory is based on an asymmetric treatment of the different 
components of the electrolyte solution. The weak coupling Debye-Hiickel approach is used in order 
to describe the monovalent salt ions while a strong coupling approach is used to tackle the polyvalent 
counterions. This combined weak-strong coupling approach effectively leads to dressed interactions 
between polyvalent counterions and thus directly affects the correlation attraction mediated by 
polyvalent counterions between like-charged objects. The general theory is specifically applied to 
a system composed of two uniformly charged plane-parallel surfaces in the presence of salt and 
polyvalent counterions. In the strong coupling limit for polyvalent counterions the comparison with 
Monte-Carlo simulations shows good agreement for large enough values of the electrostatic coupling 
parameter. We delineate two limiting laws that in fact encompass all the Monte-Carlo data. 



I. INTRODUCTION 

The properties of biological matter, soft matter and 
colloidal systems such as polyelectrolytes, charged mem- 
branes and charged micelles are in many important re- 
spects governed by the underlying properties of electro- 
static interactions In these aqueous systems the 
charged macromolecular (macroion) surfaces are always 
surrounded by neutralizing counterions as well as salt 
ions from the electrolyte solution. In simple solutions 
of monovalent salt (e.g., physiological solution of NaCl), 
the system is well described by the traditional Poisson- 
Boltzmann (PB) approach 0]. Being mean-field based, 
the PB theory is applicable only for sufficiently small 
macroion surface charge densities, low ion valencies, high 
medium dielectric constant and/or high temperatures. 
The limitations of this approach become practically im- 
portant in highly charged systems when ion-ion correla- 
tion effects start to govern the electrostatic properties of 
the system I n some sense the accuracy of the PB 
approach can be systematically improved by perturbative 
corrections in the ionic fluctuations and correlations 

along the lines of the standard approach used in the 
mean- field context. But the validity of this approach is 
severely limited since theperturbative (loop) expansion 
is weakly convergent d, 0, !H and higher-order correc- 
tions beyond the first-loop correction (Gaussian fluctu- 
ations) are very complicated [Til [l7j and in many cases 
intractable to evaluate. Thus, in many cases such pertur- 
bative corrections offer only an insignificant improvement 
over a PB approximation |3|, 0] and can therefore not be 
used to predict phenomena such as like-charge attraction 
that differ qualitatively from the mean-field predictions 

An alternative approach has been pioneered by several 
workers [fUl] using various analytical methods that allow 
for non-perturbativc treatment of correlations in highly 



charged systems. A deeper insight into the electrostatic 
properties of Coulomb fluids led to the conclusion that 
their behavior lies between two different limiting cases: 
the weak coupling (WC) limit, governed by the PB the- 
ory, and the strong coupling (SC) limit-governed by the 
so-called strong coupling theory 0, 0, Hf- In fact, these 
two limits can be established [8| as two exact limiting 
laws from a single general theory and thus allow for an 
exact treatment of charged systems at two disjoint limit- 
ing conditions. The parameter space in between can be 
analyzed by approximate methods (l9l - [2l| that cross over 
between the two limits and is mostly accessible solely via 
computer simulations @, i, |, M, H H Ex- 
act solutions for the whole range of coupling parameters 
are available only in one dimension |27fl . The WC-SC 
paradi gm has been tested extensively!!, |3, d, [ljj Il9r - 
Y2A\ l26l |27|| and was found to describe computer simula- 
tions quantitatively correctly in the respective regimes of 
validity, thus providing a unifying view of the behavior 
of Coulomb fluids. 

However, the SC theory was so far designed exclu- 
sively to treat counterion-only systems, i.e., a Coulomb 
fluid composed of only counterions in the absence of any 
salt ions |8( . Though an approximation of this type can 
be used to describe situations where a large amount of 
polyvalent counterions dominate the electrostatic behav- 
ior, it is in general unrealistic and has to be amended in 
order to deal with the complexity of real systems that 
always contain variable amounts of simple salt [28j. A 
common situation would be a system composed of fixed 
surface charges with polyvalent counterions bathed in a 
solution of simple monovalent salt [29|, This leads 
to a difficult problem of asymmetric aqueous electrolytes 
where different components of the Coulomb fluid are dif- 
ferently coupled to local electrostatic fields [3l|. Some 
components, such as polyvalent counterions, arc coupled 
strongly, whereas some, such as monovalent salt ions, 
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are coupled weakly. In systems of this type, no single 
approximation scheme that would treat all the charged 
components on the same level would be expected to work. 
Whereas the SC framework would certainly work for the 
polyvalent counterions, it would fail for the monovalent 
salt. The converse is true for the WC framework. In 
providing a consistent theoretical framework to describe 
highly asymmetric electrolytes, one is thus faced with 
a conundrum since no single approximation scheme ap- 
pears to be valid in any range of coupling parameters. 
In what follows we will build a theoretical framework 
that will allow us to selectively use different approxima- 
tion schemes for different components of the asymmetric 
Coulomb fluid. This combined WC-SC approach appears 
to bring forth all the salient features of these asymmet- 
ric systems at high electrostatic couplings that can be 
observed in Monte-Carlo (MC) simulations. 

Both the WC and the SC limits of a system composed 
of fixed surface charges with intervening mobile coun- 
terions can be obtained based on a functional integral 
or field-theoretic representation [l2j of the grand canon- 
ical partition function. The functional integral represen- 
tation of the Coulomb fluid partition function was pi- 
oneered by Edwards and Lenard in the 1960s [32[ and 
has found many applications in the modern theory of 
Coulomb fluids. It was shown by Netz Q that both 
limits can be obtained by systematic expansion with re- 
spect to a single electrostatic coupling parameter, which 
is connected closely to the plasma parameter and quan- 
tifies the relative importance of interactions between the 
mobile charges versus their interaction with a fixed ex- 
ternal charge distribution of opposite sign. It can be 
introduced in the following way: the typical distance at 
which two unit charges eo interact with thermal energy 
IzbT in a medium of dielectric constant e is known as the 
Bjerrum length, £b = eg/ (Aireeok-BT) (in water at room 
temperature, the Bjerrum length is £b ~ 0.7 nm). In 
that respect the Bjerrum length for two counterions with 
valencies +q is then q 2 Ib- Similarly, the typical distance 
at which a counterion interacts with a charged surface 
(with surface charge density of —a) with an energy equal 
to k^T is given by the Gouy- Chapman length, defined as 
Mgc = eo/(27r<7£BC r )- A competition between counterion- 
counterion and counterion-surface interactions can be 
quantitatively measured by the ratio of these two char- 
acteristic length scales, i.e. 

S = <Z 2 Vmgc = 2nq 3 £ B a/e , (1) 

which is known as the (Netz-Moreira) electrostatic cou- 
pling parameter Q. 

The WC regime is characterized by a small coupling 
parameter 5 <C 1. It is governed by the mean-field 
PB theory (obtained in the exact limit of S — ► 0) and 
weak Gaussian fluctuations around the PB solution. This 
regime is adequate for low valency counterions and/or 
weakly charged surfaces and represents the physical sit- 
uation where the width of the counterion layer /xqc close 
to the oppositely charged surface is much larger than the 



separation between two neighboring counterions in the 
layer. Thus the counterion layer behaves essentially as 
a three-dimensional gas and a collective mean-field ap- 
proach of the PB type is completely justified. For low 
electrostatic potentials the PB equation may be further 
simplified into the linear form known as the Dcbyc-Hiickcl 
(DH) equation. 

On the other hand, in the SC regime S 3> 1 (ap- 
propriate for high valency counterions and/or highly 
charged surfaces), the mean distance between counteri- 
ons, a± ~ y/qeoja (as set by the local electroneutrality 
condition) , is much larger than the counterion layer width 
(i.e., Ox/mqc ~ >/3 ^S> 1) indicating that the coun- 

terions arc highly localized laterally and form a strongly 
correlated quasi-two-dimensional layer right next to an 
oppositely charged surface. In this case, the WC ap- 
proach breaks down due to strong counterion-surface and 
counterion-counterion interactions. Since each counte- 
rion is isolated laterally in a large correlation hole and 
can move almost independently from all the others along 
the direction perpendicular to the surface, the collective 
many-body effects diminish and on the leading order, in 
the limit 5 — > oo, the behavior of the system can be 
described exactly by a single-particle SC theory which 
includes only the contribution from counterion-surface 
interactions Q. 

The range of validity of the aforementioned WC and 
SC theories (when applied to systems with finite values 
of the couplingparameter) has been explored thoroughly 
in Refs. 0, Interestingly, although both these the- 
ories are obtained as limiting results, they can still be 
applied to charged systems in a wide range of realistic 
system parameters as confirmed by extensive numerical 
simulations |, 1, 1, 0J, M MM M H3 ■ 

The above framework was introduced originally for the 
case where the mobile charges belong to a single species of 
counterions @ . In what follows we will attempt to formu- 
late a systematic description for an asymmetric system 
containing polyvalent counterions and monovalent added 
salt and delineate its various limiting behaviors in a uni- 
fying theoretical framework. In the SC limit, we shall use 
the standard SC procedure which is based on a virial and 
1 /E expansion of the partition function [1, H|, [|| , whereas 
in the WC limit, we shall perform a saddle point evalua- 
tion of the functional integral representation of the par- 
tition function p| fl2| . The main approximation that we 
introduce in the SC limit is to treat the salt ions as weakly 
coupled to all the charges in the system. Their presence 
in the system thus has a simple consequence, namely, 
it renormalizcs the Coulomb kernel into a Debye-Hiickcl 
kernel. The counterions can be treated, however, either 
as weakly or strongly coupled depending on, e.g., their 
charge valency. As we show, this approximation makes 
sense only when one deals with a highly asymmetric sys- 
tem such as the one considered here. 

In the developments as detailed in what follows, the 
WC and the SC approaches do not carry an equal weight. 
In fact we include the WC approximation in our discus- 
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FIG. 1: Schematic view of two negatively charged surfaces 
with uniform surface charge density —a and polyvalent coun- 
terions in between. Salt ions are present everywhere is space. 



sion only for reasons of completeness since it reduces to 
well-known and established results in the framework of 
the PB theory. The main thrust of this paper is thus 
devoted almost exclusively to the SC limit in the case of 
highly asymmetric systems that have thus far eluded a 
consistent statistical mechanical formulation. 



II. MODEL 

We set up a simple model system that would exhibit 
the most important features of a typical asymmetric 
Coulomb fluid. It is composed of a fixed charge distri- 
bution (macroions), /3o( r ); and polyvalent counterions of 
charge valency +q in a bathing salt solution. For simplic- 
ity we consider a simple symmetric salt with monovalent 
co- and counterions (the results can be generalized im- 
mediately to any number of salt species as long as salt 
ions stay weakly coupled to other charge species as will 
be made clear later). We shall later on apply our theory 
to the specific case where fixed charges are present on 
two plane-parallel surfaces (of infinite area S) located at 
positions z = ±a, each carrying a uniform surface charge 
density of —a, i.e., in this case 



po(r) = —a[8(z - a) + S(z + a)]. 



(2) 



The polyvalent counterions are assumed to be confined 
to the region in between the charged surfaces, whereas 
the salt ions are assumed to be present in all regions in 
space (Fig. [TJ . In what follows we may refer to the q va- 
lency (polyvalent) counterions simply as "counterions"; 
they should not be confused with salt cations which we 
shall refer to as "salt counterions" or "cations" . We shall 
neglect the effects of image charges arising from dielec- 
tric discontinuity at the bounding surfaces. The effects 
due to an inhomogeneous distribution of salt in the sys- 
tem as well as the image charges can be included in a 
straightforward manner within our formalism (26l . [33l |34| 
and will be considered elsewhere. 

So far we have invoked no constraints on the num- 
ber of polyvalent counterions, iV, between the surfaces. 



We will consider two different ensembles for counterions, 
namely, the canonical ensemble, where the number of 
polyvalent counterions is fixed in the system, and the 
grand-canonical ensemble, where polyvalent counterions 
can be exchanged between the system and the bulk ionic 
"reservoir" and N is thus determined by assuming chem- 
ical equilibrium with the bulk. Note that the salt ions 
are always treated within the grand-canonical ensemble, 
i.e., they are assumed to be in equilibrium with a bulk 
reservoir. We also introduce a dimensionless counterion 
fraction parameter 77 as the amount of polyvalent counte- 
rions between the surfaces relative to the surface charge, 
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The case 77 = represents a system with salt only, and 
77 = 1 is the case when the total charge due to counterions 
exactly compensates the surface charge. Note that in the 
canonical case the parameter 77 is fixed, whereas it is a 
function of inter-surface distance in the grand-canonical 
case. In a system with counterions only the global elec- 
troneutrality condition stipulates that r) = 1. This no 
longer needs to be the case and r\ can take any non- 
negative value when salt ions are present, which-due to 
their screening properties-turn the long range Coulomb 
potential into a short range DH potential (see below) 
and can thus take care of the electroneutrality condition 
themselves. 



III. GENERAL FORMALISM 

The total electrostatic interaction energy of a given 
configuration of the system can be written as 



W= \\ p(rMr,r>(r')drdr', 



where w(r,r') is the Coulomb kernel given by 

1 



u(r,r') 



47T££o|r — r'l ' 



(4) 



(5) 



and p(r) is the total charge density comprising both the 
fixed charge distribution as well as the charge distribution 
of mobile charges, 

i 

+ e S(B+- r) - ]T e S(Rr- r ). ( 6 ) 

i i 

Here R£, Rf and R~ are the positions of the polyvalent 
counterions, monovalent cations (salt counterions) and 
monovalent anions (salt coions), respectively. We then 
follow the standard procedure by introducing a fluctu- 
ating local potential, c/>, via the Hubbard-Stratonovich 
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transformation, which leads to the following exact func- 
tional integral representation for the grand-canonical 
partition function [1, [H| 

Z G = [vte-WM, (7) 

where j3 = l/(fcBT) and the field-functional Hamiltonian 
reads 

H[4>] = \ J J ^(r)u- 1 (r,r / )^(r , )drdr'+i J p o (r)0(r)dr 
~y J e- i/3eo ^ (r )fi(r)dr 

_ £± J e -^«o*«dr -tj- j e^^dr (8) 

Here ?i~ 1 (r, r') is the inverse Coulomb kernel, 

u-\t,t>) = -ee Q V 2 6(r-r>), (9) 

and A c and A± represent the fugacities of polyvalent 
counterfoils and salt ions. The geometry function f2(r) 
specifies the region accessible to polyvalent counterions; 
for the plane-parallel system in Fig. [TJ we have fi(r) = 1 
for \z\ < a and zero otherwise. The effective Hamiltonian 
([8]) involves nonlinear terms that result from integrating 
out the microscopic degrees of freedom due to counteri- 
ons and salt ions and elude an exact analytical treatment. 
We shall thus concentrate on the asymptotic behavior of 
the system in the WC and SC limit where we show that 
the problem can be handled analytically for an asymmet- 
ric system as considered in this work. 

In the WC limit (for both the counterions as well as the 
salt ions), the above functional integral is dominated by 
the saddle point of the Hamiltonian H[<fi] which leads 
to the standard PB theory In an asymmetric system 
with large q the monovalent salt ions can be treated in 
the linearized DH limit as we shall discuss in Section ITVl 
below. 

In the SC limit, the WC approximation breaks down 
and one should apply a different approach that we shall 
refer to as the SC dressed counterion theory. In this ap- 
proach, due to the high asymmetry for large q, the mono- 
valent salt ions can be treated again as weakly coupled 
on the DH level even though the counterions are treated 
as strongly coupled. This situation will be considered in 
Sections IVllVIII It will be shown that these asymptotic 
theories constitute two limiting laws that in fact bracket 
all the Monte-Carlo data. 



IV. WEAK COUPLING LIMIT: MEAN-FIELD 
PB THEORY 

As already stated in the Introduction we include the 
discussion of the WC limit only for reasons of complete- 
ness. The WC limit does not venture outside the stan- 
dard PB theory except marginally in the fact that we 



derive a limiting form of the standard PB equation valid 
in the case of a highly asymmetric ionic system when 
salt ions are monovalent but counterions have a larger va- 
lency q, adjusting other physical parameters of the theory 
(such as the surface charge density) in such a way that 
the coupling parameter for salt ions as well as counterions 
remains by assumption small. The thus constructed PB 
theory in fact represents the formal weakly coupled limit 
of the strongly coupled theory of dressed counterions to 
be introduced later. 

As is well known from the previous work on the 
counterion-only systems, the WC limit corresponds to 
the saddle-point evaluation of the partition function in 
the limit of a vanishing coupling parameter @, 0, Hi, i-e. 
via the condition 

5H[<f>]/6<f>U SP = 0. (10) 

This subsequently leads to the standard PB equation de- 
termining the real-valued mean-field electrostatic poten- 
tial ip = iifisp, that is 

- ee V 2 V = e qAMr)e-P eoq,f ' + 

+ eoA+e-^ - e A_e M = 
= e qn c (r) + e n+(r) - e n_(r). (11) 

Note that the terms on the right hand side represent 
the PB charge densities due to valency q counterions, 
s-oQ n c {r) , and the monovalent salt ions, ±eon±(r), re- 
spectively. 

The boundary conditions at both surfaces require the 
continuity of the potential ip and jump in its first deriva- 
tive due to the surface charge, 

^'(±a)i„-V'(±a)out = T — ■ (12) 

ee 

There are two ways to proceed now depending on the 
thermodynamic ensemble chosen for the counterions. As 
noted before, the salt ions are assumed to be exchange- 
able with the bulk reservoir and are thus treated always 
within a grand-canonical description. But the number of 
counterions may be taken to be fixed or determined via 
an equilibrium with the bulk. 

A. Canonical ensemble 

In the canonical ensemble the counterions are not in 
equilibrium with a bulk reservoir. Their number in the 
slit is fixed and given by 

A c J n(r)e- peoq *dr = N. (13) 

As already stated, we will express the number of coun- 
terions through the parameter 77, Eq. ([3]). On the other 
hand, the salt ions are in equilibrium with a bulk reser- 
voir that contains equal amounts of both positive and 
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negative ions, which implies that they have equal fugac- 
itics A + = A_ = ?iq, giving the DH screening parameter 
k (inverse "screening length") as 

k 2 = 8ir£ B n . (14) 

At this stage it is convenient to introduce dimensionlcss 
quantities, namely 

f = r//i G c, K = Kfi GC , u = (3e qip, (15) 

where /xgc is the Gouy-Chapman length as defined before 
and u is the dimensionless potential. 

The dimensionless PB equation now assumes the form 

u" = K 2 qsinh- +Ce~ u . (16) 

q 
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FIG. 2: (Color online) Rescaled canonical PB pressure (|18p 
as a function of the rescaled inter-surface half-distance a for 
97 = 1 and different values of the screening parameter k as 
shown on the graph. 



In the case of a highly asymmetric ionic system with q ^> 
1, one can linearize the salt ions term in the PB equation, 
leaving the counterion term in the non-linear form. This 
assumption is valid if the dimensionless potential u/q is 
small enough, leading to the PB equation in the space 
between the surfaces \z\ < a of the form 



u" = k 2 u — Cc 



(17) 



which is independent of the valency q. The constant C 
can be evaluated when one stipulates the fixed amount of 
counterfoils, eq. (|13p . which gives C f° a exp(— u) dz = 4?/. 
In the same limit the PB equation outside the surfaces, 
\z\ > a, has the form u" = k 2 u, which yields the well- 
known solutions u{z) = uo exp(±Kz). 

The interaction pressure, p, between the bounding 
surfaces is given by the difference of the ion concen- 
trations at the mid-plane (2 = 0), where the mean 
electric field vanishes, and the bulk concentration, i.e., 
pp = n+(0) + n_(0) + n c (0) — 2n , which leads to the 
dimensionless expression 



1 ~i 1 

p = —K U 



(0) + ^Cc^°\ 



where we have defined the dimensionless pressure 
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(18) 



(19) 



As evident from the above equation, the pressure can 
never be negative, i.e., the effective interaction between 
two like-charged surfaces is always repulsive within this 
type of mean-field approach (this statement is true also 
within the full nonlinear theory) (3f| . The canonical PB 
equation (|17|) can be solved numerically for the dimen- 
sionless potential u. The results for the PB pressure as 
a function of the inter-surface separation are shown in 
Fig. [5] for r} = 1 and several different values of k. The 
pressure depends only weakly on the screening parameter 
due to dominant osmotic contributions from counterions. 



B. Grand-canonical ensemble 

In the grand canonical ensemble the counterions too 
are in equilibrium with the bulk. In what follows the 
bulk concentration of counterions cq will be expressed by 
the parameter 



X 2 = 8irq 2 e B c - 



The amount of positive n Q y 1 and negative 



(20) 



<■(-) 



salt 

ions in the bulk containing also counterions is not equiv- 
alent. Here the electroneutrality condition demands that 

*(+) , *(-) 
?i v + qco = ?i v . 

We allow bulk counterions to enter only the space be- 
tween the surfaces (i.e., inside the slit \z\ < a), but they 
are not allowed to pass through the bounding surfaces 
and enter the "outside" region \z\ > a. In contrast, salt 
ions are in equilibrium both with the bulk as well as the 
region \z\ > a. In this situation the Donnan equilibrium 
is established between the bulk reservoir and the "out- 
side" region, giving therefore 



*(±) 



n e 



(21) 



where no is the salt concentration infinitely far from the 
charged surfaces in the "outside" region and tpu is the 
corresponding Donnan potential (i.e., the potential dif- 
ference between the bulk and the region \z\ > a where 
the potential is in fact assumed to be zero infinitely far 
from the charged surfaces). From the above equations 
we obtain 



. , u D qc 
smn — = 



X 



2uq 2n 2 q 



(22) 



for the dimensionless Donnan potential ud = fteoqipD, 
which, to the first order (appropriate for the linearized 
theory), gives «d = x 2 /2k 2 . Note that we have defined 
the screening parameter k in terms of the concentrations 
in the "outside" region as k 2 = 87t£b^Oi which is identical 
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FIG. 3: (Color online) Rescaled grand-canonical PB pressure 
(|25|l as a function of the rescaled inter-surface half-distance a 
for k = 0.5 and different values of the parameter \ as shown 
on the graph. 



to the one used in the canonical case. We also have A + = 
A_ = n and A c = coe Un . 

The grand-canonical PB equation in the slit region 
takes the following form for the dimensionless potential, 



K 2 g sinh--ix 2 e- {u ~ UD) . 
q 



(23) 



After linearization of the salt part (first term on the right 
hand side of Eq. (|2"3")l ) in the case of a highly asymmetric 
system with large q, the PB equation assumes the form 



u" = K 2 U ■ 



h 2 



-(u—ud) 



(24) 



In the region \z\ > a, the PB equation has the same 
form and solutions as in the canonical case, see Section 
IIV Al We furthermore stipulate the same boundary con- 
ditions as in the canonical case. The inter-surface pres- 
sure is again given by the difference of the ion concentra- 
tions at the mid-plane (z = 0) and in the bulk so that 

f3p = n+(0) + n_(0) + ?i c (0) - n* Q (+) - n*^ - c . The 
dimensionless pressure then assumes the form 

p = h 2 (u 2 (0) - u 2 D ) + l£ 2 [e-W°)--) - 1], (25) 
where we have introduced the dimensionless parameter 



(26) 



One can again show that the grand-canonical pressure 
is always positive within the PB theory, just as in the 
canonical case [35|. The PB equation in the grand- 
canonical ensemble, Eq. (|2"4")l , is solved numerically. The 
results for the PB pressure are shown in Fig.[3]for k = 0.5 
and several different values of x- Note that the grand- 
canonical pressure remains finite also in the limit of van- 
ishing inter-surface separation as in this case all the coun- 
tcrions can be expelled from the inter-surface region. 



V. SC DRESSED COUNTERION THEORY: 
FORMALISM 

We now turn our attention to the main part of this 
paper, namely the case with strongly coupled counteri- 
ons. Proceeding from the Hamiltonian ©, we introduce 
the following approximation: we assume that only the 
polyvalent countcrions are strongly coupled, whereas the 
monovalent salt ions are weakly coupled to other charges 
in the system. Since salt ions are in equilibrium with 
a bulk reservoir, the fugacities of salt counterions and 
coions are the same, A + = A_ = no, and have exactly 
the same meaning as in the previous Section. 

The salt terms (the last two terms) in Eq. (|5|) can 
be expressed as cos f3eo<fi(r) and for a highly asymmetric 
system they can be expanded up to the quadratic order 
in the fluctuating potential. Thus up to an irrelevant 
constant we remain with 

H eS [<j>} = \ JJ 0(i> D ^(r,r')0(r')drdr' + i / p o (r)0(r)dr 

_ ^ [ c- i/3eo ^ r }ft(r)dr. (27) 
P J 

This procedure therefore yields an effective Hamiltonian 
for a "counterion-only" system but with the proviso that 
the inverse Coulomb kernel is replaced by the standard 
inverse Dcbyc-Huckcl kernel 



DH 



(28) 



where k 2 = ^■kI-qtl^. We have thus effectively integrated 
out the salt degrees of freedom leading to a renormalized 
interaction potential between all the remaining charge 
species of the screened DH form, 



u DH (r,r') = 



e -re|r-r'| 
47T££o|r — r' 



(29) 



One can thus drop any reference to explicit salt ions and 
infer thermodynamic properties of the original system by 
analyzing it as a system composed of only "dressed coun- 
terions" (as well as fixed external charges) that interact 
via a screened DH pair potential. We term this approxi- 
mation scheme in the SC limit as the SC dressed counte- 
rion theory. Our SC analysis thus proceeds in the same 
way as for the counterion-only systems @ except that the 
interactions between the charges are dressed (Sections IVII 
and Em . 

We note again that on this level the salt ions are al- 
ways considered to be weakly coupled to all the other 
charges as well as between themselves. Any Bjerrum 
pairing [H, H3 or even electrostatic collapse of the salt 
or formation of salt-counterion complexes [38( is thus be- 
yond the approach developed here. The regime of param- 
eters where these effects may become important and the 
regime of validity of the SC dressed countcrion theory 
will be determined later. 
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VI. SC DRESSED COUNTERION THEORY: 
CANONICAL ENSEMBLE 

We first discuss the SC limit within a canonical en- 
semble for the polyvalent countcrions. This means that 
the number of counterions in the system is fixed. The 
analysis in this case is very similar to the traditional SC 
approach Q. In the SC limit, we expand the grand- 
canonical partition function associated with the dressed 
counterion approximation, Eq. (|27p. to the first or- 
der in counterion fugacity, A c , and then perform an 
inverse Legendre transformation (where the fugacity is 
mapped to the number of countcrions, N, via the rela- 
tion A c dlnZ(j /dA c = N) in order to obtain the canonical 
SC free energy [H, 0, IH, HiJ ■ For the system under con- 
sideration we find 



where we have introduced 



j3F N = (3W o -Nhx / e- pWo '^dr 



(30) 



where the first term is the screened interaction energy of 
fixed charges (macroions) 



Woo = r 



drdr'/9o(r)MDH(r, r')p (r'), (31) 



and the term in the exponent is the single-particle in- 
teraction energy of the dressed counterions with fixed 
macroion charges 



W 0c (r) = e q / drVo(r') u DH(r, r ') 



(32) 



The SC attraction between like-charged macroions stems 
from the second term in Eg. (1301) . which contains the 
counterion-induced effects |3S,l33|- 

We now focus on the case of two charged planar sur- 
faces as defined in Section [TTJ In the canonical ensemble, 
we stipulate that the number of countcrions N between 
the bounding surfaces remains fixed even as the distance 
between the surfaces is changed. The number of countcri- 
ons, N, may be arbitrary and, unlike the counterion-only 
case Q , its value is not determined by the electroneutral- 
ity condition. The DH salt takes care in neutralizing all 
charges as follows from the observation that the dressed 
interaction potentials are short ranged. 

Expressing the interaction energies with dimensionless 
quantities, Eq. (|15j) . and defining the rescaled energy 

Woo/S = (2££ Woo)/(c 2 /xgc<S') per rescaled area S = 
S/[Aq C , we obtain 



Woo/S=\e- 2ki , 

K 

j3W 0c (z) = ^e~ kii cosh kz, 



(33) 
(34) 



and thus the following expression for the SC free energy 
of the plane-parallel system, 



Pn/S 



1 

— e 

K 



2ka 



2r)hil(a), 



(35) 



1(a) 



( 2 - 
cxp — e 



cosh kz ] dz 



(36) 



The first term in the free energy corresponds to the usual 
salt-mediated DH repulsive interaction between the sur- 
faces, and the second one is the contribution of coun- 
tcrions and is thus proportional to their amount in the 
slit expressed by the dimensionless parameter rj, Eq. ([3]). 
The dressed counterion (number) density profile is ob- 
tained in the form of the Boltzmann factor of the single- 
particle potential Wq c (t) (i.e., the integrand in Eq. ([36)) ). 
and follows as 



/ 2 

c(z) — C cxp( — e" 



' cosh kz 



(37) 



The normalization constant C can be deduced from the 
known amount of counterions, rj. The canonical SC the- 
ory contains two free physical parameters, namely, the 
salt screening parameter k and the parameter rj. The 
model itself additionally contains the coupling parame- 
ter, S, which is considered as S — > 00 within the SC 
theory. We will use MC simulations to compare the re- 
sults of the SC theory with those obtained numerically 
at finite values of S. 

The dimensionless pressure acting on each surface 
can be obtained from the free energy via the stan- 
dard thermodynamic equalities in the form p = 
— \ J '§)/ '(da), thus yielding 



P 



I' (a) 
1(a)' 



(38) 



where the prime denotes the derivative with respect to 
the argument. The dependence of the SC pressure on the 
inter-surface separation and the screening length will be 
discussed in detail later (Fig. [Sk) ■ Here we briefly con- 
sider its limiting behavior at small and large separations. 

In the limit of small separations ko < 1, (i-e., small 
compared to the screening length k" 1 ), the dimensionless 
pressure reduces to 



p~? + (1-2*7). 
a 



(39) 



The leading order term corresponds to the ideal-gas os- 
motic pressure of counterions which is repulsive and dom- 
inates over other electrostatic contributions. In the limit 
of large separations tea > 1 (i-e., large compared to the 
screening length k -1 ), the pressure behaves as 



~ ~ 1 - v 
a (ka) 2 



(40) 



which indicates that at large separations the counterions 
also behave as an ideal gas since all their electrostatic 
interactions arc effectively screened out at several n~ l 
distances. Hence, only the repulsive osmotic contribution 
remains. 
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TABLE I: Illustrative examples of possible physical parame- 
ters corresponding to a given set of values for the electrostatic 
coupling parameter E (in water at room temperature with 
Ib ~ 0.7 nm). no"™ shows the minimal salt concentration 
in order to satisfy the criterion (|42[) . The salt concentration 
no corresponds to the inverse dimensionless screening length 
ii — 1 and the bulk polyvalent counterion concentration Co 
corresponds to x = 1. For 3 = 500 water as the dielectric 
medium (e = 80) is replaced by ethanol (e = 30). 



The canonical dressed counterion theory can be used 
to derive the well-known SC pressure in the no-salt limit 
po = 1/a— 10) where only (polyvalent) counterions are 
present without any salt ions. This limit is recovered by 
letting k — > and setting rj — 1 in order to satisfy the 
electroneutrality condition. 



A. Regime of validity of the canonical SC dressed 
counterion theory 

The regime where the SC dressed counterion theory is 
applicable follows from a simple criterion that was de- 
rived and confirmed by MC simulations 0, 0, Q. The 
pivotal statement in deriving this criterion is that the 
separation between nearest-neighbor counterions must be 
much larger than the separation between the surfaces, 
a_L ^ o,. In dimensionless units the separation between 
neighboring counterions is approximately a±_ ~ •v/3/r?. 
The inter-surface separation at which our SC theory 
should be valid is then 



a < 



(41) 



This indicates that the regime of validity of the SC theory 
increases as the coupling parameter becomes larger. In- 
terestingly, the regime of validity of the theory increases 
also as r) becomes smaller (at fixed coupling). In partic- 
ular, for rj < 1 (i.e., when the amount of the bare charge 
due counterions is less than the bare fixed charge on the 
macroions) the SC dressed counterion theory is expected 
to hold in a wider range of separations as compared with 
the original no-salt SC theory Q. 

On the other hand, at very large inter-surface separa- 
tions, where most of the electrostatics is screened out, 
the interaction between counterions becomes negligible 
and the SC theory of the dressed counterions, Eq. (f55l) . 
retains its validity again. We therefore expect that the 
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FIG. 4: (Color online) a) Normalized SC dressed counterion 
density profile between two charged planar surfaces within the 
canonical description, Eq. (|37|l . for half-distance a — 1 and 
various screening parameters k. The normalized SC density 
is independent of the parameter n. b) The theoretical SC 
density profile (solid line) is compared with MC simulation 
data (symbols) for R = 1, n = 1, S = 4 and two different 
coupling parameters E as used in the simulations. 



SC for dressed counterions ceases to be valid only at in- 
termediate separations. 

The validity of the DH-type linearization that we have 
used to derive the dressed interaction potentials is also 
limited by stipulating that the dimensionless DH poten- 
tial itself is always small enough. Specifically this sets 
the value of the dimensionless DH potential at the two 
boundaries to be well below unity, which translates to 
the condition e$K 3> 2n^a, or 



ii > 1/q, 



(42) 



where q is the counterions valency. Table U shows illus- 
trative examples of possible physical parameters where 
the above criterion is fulfilled. 



B. MC simulations: Canonical ensemble 

In order to assess the validity of the SC dressed counte- 
rion theory, we performed MC simulations and compared 
them with the analytical results. Simulations cover a 
wider range of coupling parameters and can thus fill the 
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FIG. 5: (Color online) Rescaled canonical pressure as a function of the rescaled inter-surface half-distance a. In (a) the SC 
results obtained from Eq. (|38|l are shown for 77 = 1 and different values of the screening parameter k. In (b) and (c) the SC 
result (solid line, Eq. (1381) ) and the PB prediction (dashed line, Eq. (JTHJ) ) are compared with MC simulation data (symbols) 
obtained for different values of the coupling parameter and at fixed k and rj as indicated on the graphs. 




gap between the WC and the SC coupling limits investi- 
gated analytically in the preceding Sections. The simu- 
lation details are described in Appendix [5] 

First we compare the counterion density profiles ob- 
tained from the simulations with the analytical result, 
Eq. (|57|) . As seen in Fig. HI the counterion density pro- 
file is a convex function of the inter-surface position so 
that the counterions are preferentially gathered near both 
charged surfaces due to electrostatic attraction. For large 
screening or large separations, na 3> 1, the typical coun- 
terion density decay length from each surface is In 
contrast, at low screening or small separations, na <C 1, 
the typical decay length is \J «~ Vgc- Note that the den- 
sity profile has non-monotonic variations as the screening 
is increased. For vanishing k Q as well as for very large 
k the SC dressed counterion density profile becomes ho- 
mogeneous, p(z) = const. The agreement with MC data 
is excellent already at 3 = 50 as shown in Fig. 2J 

Next we consider the inter-surface pressure. The ana- 
lytical pressure-distance plots from Eq. (pt5| are shown 
in Fig. [5ji for a few different screening parameters. At 
small separations the pressure is always repulsive due to 
the large osmotic component of counterions stacked in 
the small region between the surfaces and is indepen- 
dent of the screening, Eq. (|59")l . By increasing the inter- 



surface separation, the osmotic pressure as well as the 
electrostatic counterion-induced attraction (embedded in 
the second term of Eq. (|30|) or ([38]) ) are reduced. The 
repulsive osmotic pressure decays approximately as 1/a 
for small distances, whereas attractive electrostatic con- 
tributions decay with the length scale n~ l . Thus, when 
k is small enough, the attractive electrostatic component 
may remain large enough at intermediate distances to 
overcome osmotic repulsion and hence a net attraction 
can appear. 

As clearly seen form Fig. [5^,, the larger the screening 
parameter k the smaller will be the attraction at inter- 
mediate separations. For large enough k the attraction 
disappears altogether and the pressure reverts to its lim- 
iting form, Eq. (|40[) . For large inter-surface separations 
the pressure thus always ends in a repulsive regime due 
to slow (1/a) decay of the osmotic contributions as com- 
pared to rapid (e~ Ka ) decay of the electrostatic contribu- 
tions. 

Comparison with MC results in Figs. [5}d and c shows 
expectedly that the agreement becomes better as the cou- 
pling parameter 3 becomes larger. The agreement is also 
better for a smaller fraction of counterions, 77, in the slit 
(cf. Figs. [5] and which also follows from Eq. (jUJ). At 
large separations or strong screening the electrostatics is 
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FIG. 7: (Color online) Van der Waals type iso-inverse screen- 
ing length curves are shown along with the corresponding 
Maxwell construction for the interaction pressure between 
charged plates within the canonical description for dressed 
counterions. Here r) = 1 and the dimensionless inverse screen- 
ing length varies in the range from k — 0.4 to k — 0.56 in in- 
tervals of 0.027 (from bottom to top). Red line represents the 
critical curve. The critical point is located at approximately 
p c = 0.092, a c = 3.07 and k c = 0.546. 



screened out and the SC dressed counterion approxima- 
tion becomes valid again. 

The validity of the SC approximation at large separa- 
tion or strong screening is a consequence of the dressed 
counterion theory and is not obtained in the standard 
SC theory with counterions only, which remains valid 
only in the regime of small separations 0, 0, 13 • I n the 
dressed counterion approach at high screening or large 
separations the counterions behave as "ideal" particles 
and do not interact anymore with fixed charges in the sys- 
tem leading to a repulsive ideal-gas pressure regime with 
the 1/a signature dependence on the inter-surface separa- 
tion. This dependence is not obtained in the counterion- 
only SC theory where the interactions remain attractive 
and long range at large separations. We may thus con- 
clude that the SC dressed counterion theory captures the 
physics both at large and small separations and requires 
improvements only at intermediate distances. 

The interaction pressure in the canonical system thus 
always possesses repulsive branches at small and large 
separations and can show non-monotonic behavior in be- 
tween as seen in the figures. In fact, for certain values 
of the parameters the interaction pressure shows a van 
der Waals-like loop which could suggest a coexistence 
regime between two different "phases". The difference 
between the standard van der Waals isotherms and the 
case encountered here is that the curves are not isotherms 
but are rather iso-inverse screening length or iso-ionic 
strength curves. From a thermodynamic perspective, this 
would indicate a dense phase identified with a small inter- 
surface separation at equilibrium with an expanded phase 
with a large separation. Such van der Waals-type loops in 
the interaction pressure between interacting charged sur- 
faces have been seen in other contexts before [40|. The 
phase coexistence can be demonstrated by means of a 



Maxwell construction analysis as shown in Fig. [7] In an 
experiment that can only probe stable equilibrium states 
in the system, as is the case in osmotic stress experi- 
ments and in general also in the surface force experiments 
[2&| , one would measure exactly these type of interaction 
pressure versus separation curves that are in agreement 
with the appropriate Maxwell construction. The binodal 
or the coexistence curve, that delimits the region in the 
pressure-separation dependence where Maxwell construc- 
tion is feasible, ends at a critical point corresponding to 
the largest amount of salt in the system that still leads 
to a non-monotonic dependence of the pressure on the 
inter-surface separation. Above this limiting salt concen- 
tration the interaction pressure is purely repulsive. For 
the case with 77 = 1 we numerically obtain the critical 
point as p c = 0.092, a c = 3.07 and k c = 0.546 in rcscalcd 
units. 

It is interesting to note that this type of interaction 
pressure equilibria corresponding to abrupt transitions 
from one equilibrium separation to a different one are in- 
deed observed between strongly charged macromolecules 
in the presence of polyvalent counterions and monovalent 
salt. Most notably osmotic stress experiments on DNA in 
the presence of trivalent CoHex counterions and 0.25M 
NaCl salt show that for intermediate Co 34 " concentra- 
tions (8 mM and 12 mM) there are abrupt transitions 
from one repulsive osmotic pressure branch to another 
one (29|. Similar features are discerned even for a diva- 
lent counterion Mn 2+ at various concentrations or tem- 
peratures (30j . 

VII. SC DRESSED COUNTERION THEORY: 
GRAND-CANONICAL ENSEMBLE 

So far we have considered only the canonical ensem- 
ble where the system contains a fixed amount of dressed 
(polyvalent) counterions. In this Section we focus on the 
grand-canonical ensemble for the SC dressed counterion 
theory. In the grand-canonical ensemble, the system is 
in equilibrium with an external reservoir of dressed coun- 
terions with the bulk concentration cq. Depending on a 
chemical potential (or the fugacity A c ), the system can 
adjust the number of counterions by exchanging them 
with the reservoir. The number of counterions, N, in the 
inter-surface region thus varies with the distance between 
the surfaces. 

Here we will derive the grand-canonical free energy by 
using the canonical free energy, Eq. (|35[) , discussed in the 
previous Section. The grand canonical partition function 
is defined as 

N=0 

where Zn = c~ ,3:Fjv is the canonical partition function. 
The grand-canonical free energy is then obtained as 

/?$ = (3W Q0 - AJ(a), (44) 
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where Woo is the screened interaction energy of the 
charged surfaces (in the absence of counterions), Eqs. 
(|3Tj) and (J33J), and 1(a) is given by Eq. ([36]). We can 
evaluate the average number of counterions in the slit as 



w =-^L =a ' j(o) ' 



(45) 



and the pressure is given in terms of this latter quantity 

as 



Psub = 



1 



25 



— Nk B T— \nl(a) 

oa oa 



(46) 
(47) 



In order to satisfy the chemical potential equivalence 
with the reservoir, counterions must be distributed in the 
inter-surface region according to the Boltzmann weight, 



c(r) = c c 



(48) 



Here Wo c (r) is the single-particle surface-counterion in- 
teraction energy, Eqs. ([32]) and (j34)) . The rescaled 
dressed counterion density in the inter-surface region is 
then obtained as 



c(z) = Co exp 



cosh 



(49) 



Again this single-particle form of the density distribu- 
tion is expected in the SC limit as counterion-counterion 
interaction effects are negligible on the leading order 
[E EL HI ■ The density profile in the grand canonical (Eq. 
(|4"9"| ) and canonical ensembles (Eq. (f3"T)) ) have similar 
forms, except that the normalization constants are dif- 
ferent. In the present case, Co is the bulk concentration 
of dressed counterions and is therefore a free parame- 
ter. Integrating the counterion concentration in the inter- 
surface region, we get the average number of counterions 
as 



N 



c(r) dV = S^gc c I(a). 



(50) 



The expression (|47j) for the pressure includes only the 
contribution due to the free energy of the ions in the sys- 
tem. In order to evaluate the net pressure acting on each 
surface, one must also account for the pressure exerted 
by the counterions in the bulk (reservoir), i.e. 



PciO = c k B T. 



(51) 



The net pressure, p = p su b — Pdo, is then obtained as 



p = c- 2k « + lx 2 ll'(a)~2] 



(52) 



in rescaled units, where x 2 = ^Q 2 ^b c o Mgc The g ran d- 
canonical theory therefore contains two free physical pa- 
rameters, i.e., ii and x, which are given by the bulk salt 



and polyvalent counterions concentrations, respectively, 
as well as the rescaled half-distance, a. 

We can re-express the above pressure explicitly as 



V 



3 -/3Woc(a)_ 



1- 



c -S(5+5) e -/3Wo c (i) d ~ 



(53) 

which allows for a simple physical interpretation. The 
first term is the usual DH repulsion between the plates, 
exactly the same as in the canonical case, Eq. (f3"5|). The 
remaining part, proportional to x 2 , is the dressed coun- 
terions contribution, i.e. 



Vc 



1-2 

4 X 



-^Woo(S) 



1 



-K(i+5) e -/3W 0c (i) d - 



(54) 

The first two terms in the brackets correspond to the dif- 
ference between counterion osmotic pressures inside the 
system and in the bulk, which is hence repulsive. The last 
term in the brackets contains the attractive electrostatic 
interaction mediated by dressed counterions between the 
two surfaces. 

The pressure at small separations, na <C 1, varies lin- 
early with the distance as 



p~l 



1 -X 2 (e^-l)-(2k + x 2 e^)d. 



(55) 



Note that, in contrast to the canonical case, the grand- 
canonical pressure reaches a finite value as the separa- 
tion tends zero, since in this limit all counterions are 
expelled from the inter-surface region. For large separa- 
tions, > 1, the pressure approaches zero asymptoti- 
cally from the negative side 



V 



i*Lae- 2 ™. 
2 K 



(56) 



The above limiting form involves contributions from elec- 
trostatic as well as entropic effects that both tend to zero 
at large separations. 

In contrast to the canonical case, discussed in the 
previous Section, the grand-canonical case has no 
countcrion-only analog that could be obtained simply by 
taking the limit n — > 0. The reason is that the amount of 
counterions between the surfaces in the grand-canonical 
ensemble is controlled by the inter-surface potential and 
hence by the amount of all ions between the surfaces. 
Since in the limit k — > 0, only counterions are present, the 
electroneutrality which implies a fixed amount of coun- 
terions cannot be generally achieved. The conclusion 
is that the grand-canonical ensemble cannot be used to 
study the low-salt (or zero salt) limit, where a canonical 
description is more appropriate. 



A. Regime of validity of the grand-canonical SC 
dressed counterion theory 

The amount of counterions between the surfaces is 
fixed in the canonical case and can be expressed in terms 
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of the parameter 77, Eq. (|3|). In the grand-canonical case 
counterions can be exchanged with an external reservoir. 
Defining the grand-canonical fractional amount of coun- 
terions, 77GC1 m the slit in accordance with the definition 
([3]), we find that 



»7Gc(a) = 1(a) 



(57) 



In order to obtain a validity criterion for the grand- 
canonical SC theory, we follow the same general argu- 
ment as used in the canonical case, i.e., the SC de- 
scription is expected to hold when d <C aj_, and hence 
a <C ^/S/t7gc- As discussed before, ?7gc now depends on 
the separation a. Using Eq. (f5"T|) . we obtain the following 
criteria 



5V3 

5< ^273 C ~ 2/3K forwo-Cl, 



a < 



S x/3 
W 



(58) 



for na ^ 1, 



which clearly merge into a single criterion if k 1. 
These criteria determine the range of parameters where 
the grand-canonical SC dressed countcrion theory may be 
applied at finite couplings, e.g., when the theory is com- 
pared with numerical simulations. Note that the range of 
the validity of the grand-canonical theory is more limited 
than that of the canonical theory as the dependence on 
the coupling parameter, S, turns out to scale as S 1 / 3 in 
contrast to the usual H 1 / 2 behavior found in the canoni- 
cal case or within the standard no-salt SC theory [!, 0, [|| . 
But the regime of validity of the theory increases as the 
bulk concentration of counterions is lowered, which is 
similar to the trend found within the canonical theory. 



B. MC simulations: Grand-canonical ensemble 

An interesting quantity that can be investigated in 
the grand-canonical case is the amount of counterions, 
?7gc! between the surfaces as the inter-surface separa- 
tion, 2a, is varied. Since t]qc depends linearly on the 
amount of counterions in the bulk, the quantity ?7gc Ix 2 
depends only on the screening parameter k in addition to 
the inter-surface separation, 2a. In Fig. [5J we show the 
"charging" curve, fyec/x 2 = ?7Gc/x 2 ( a )i which reflects 
some of the properties of the counterion density profile. 
The amount of counterions starts to increase linearly at 
small separations, na <C 1. This is due to the fact that at 
small separations counterion density profile is almost con- 
stant in the inter-surface region and thus the number of 
counterions is proportional to the inter-surface distance. 
On the other hand, at very large separations, na ^> 1, 
the counterion density profile is again very flat in the 
middle of the inter-surface region, which again results in 
a linear dependence of the amount of counterions on the 
inter-surface distance as shown in the figure. 
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FIG. 8: (Color online) Charging curves of the grand- 
canonical system: the rescaled amount of counterions in the 
slit, 7/gc/x 2 ! i s plotted as a function of the rescaled half- 
distance a for various screening parameters k. The theoretical 
SC values of t/gc/x 2 (lines) are independent of \. Symbols 
correspond to MC results for k = 0.5 which indicate better 
agreement with the theoretical prediction for larger H and 
smaller \- 



It is interesting to note that for screening parame- 
ters smaller than k < 0.575, the charging process is 
non-monotonic in an intermediate range of separations. 
This is a consequence of the interplay between the inter- 
surface separation and the resulting "dressed" electro- 
static potential. In the non-monotonic part of tjqc, the 
increase in the separation results in a larger drop of the 
resulting DH potential between both surfaces. The corre- 
sponding potential drop then results in a larger decrease 
in the number of counterions than would follow from the 
increase of the separation itself and hence counterions 
are expelled from the inter-surface space as a result of 
the increase in the inter-surface distance! 

As seen in Fig. [51 the theoretically predicted value of 
7?gc gives the upper bound limit and the MC results at 
finite coupling parameter are always smaller. The agree- 
ment between the theory and simulations is better for 
larger coupling parameter S. In fact as S is decreased 
the counterion-counterion repulsions become more dom- 
inant and the amount of counterions in the slit is de- 
creased. The agreement is also better for smaller concen- 
tration of counterions (smaller x), which implies larger 
nearest-neighbor distances between the counterions and 
thus smaller counterion-counterion repulsions. These ob- 
servations agree qualitatively with the criteria (|5B)>. 

There are some qualitative differences between the SC 
dressed counterion pressure in the two ensembles. Most 
notably, the grand-canonical pressure (shown in Fig. (5^) 
does not diverge and has a finite limiting value at small 
separations, Eq. (|55p . while the canonical pressure di- 
verges due to the diverging contribution of the confine- 
ment entropy of counterions. Recall that the total pres- 
sure in the grand-canonical case, Eq. ([53]). is composed 
of a contribution due to the surface-surface repulsion, 
which is proportional to exp(— 2«a), and a contribution 
due to dressed counterions. As seen from Eq. (|5H) . the 
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FIG. 9: (Color online) Rescaled grand-canonical pressure as a function of the rescaled inter-surface half-distance a. In (a) the 
SC results obtained from Eq. (|53[) are shown for k = 1 and different values for the parameter x, which is related to the bulk 
concentration of dressed counterions as defined in the text. In (b) the rescaled counterion contribution, p c /x 2 (Eq. (|54p ). to 
the inter-surface interaction pressure is shown for different values of the screening parameter k as indicated on the graph. In 
(c) the theoretical SC (solid line, Eq. (|53[) ) and PB (dashed line, Eq. (|25[) ~) results are compared with MC results (symbols) 
obtained for different values of the coupling parameter and at fixed % = 1 and k = 0.5. 



latter contribution, p c , is proportional to the amount 
of counterions in the bulk, x 2 - Hence, p c /x 2 depends 
only on the inverse screening length k and the separa- 
tion a. The behavior of this contribution is shown in 
Fig. [2b which shows a minimum at rather small separa- 
tions whose depth decreases rapidly with increasing the 
screening parameter k. Whether the total pressure ex- 
hibits an attraction regime (with negative total pressure) 
at a given screening parameter k depends on the amount 
of counterions in the system through x 2 ■ 

In Fig. [SJ;, we compare the SC pressure with the MC 
simulations in the grand-canonical case. There is a good 
agreement at small separations in accordance with the 
criteria (|58|) . At large separations better agreement may 
be obtained by taking a larger value of the coupling pa- 
rameter. In the grand-canonical case, the interaction os- 
motic pressure decays exponentially, Eq. (|56[) , and there- 
fore electrostatic effects, which also decay exponentially, 
remain important even at large separations. 

As seen the simulation results are again bracketed by 
the two limiting analytical theories of weak and strong 
coupling within the dressed counterion scheme and thus 
agree with the general feature obtained before [H, |3, Q 
that the WC and SC limits in fact establish the upper 
and lower bounds for the interaction pressure between 
charged surfaces. 

VIII. DISCUSSION 

A. Comparison between canonical and 
grand-canonical ensembles 

In the preceding Sections, we analyzed the behavior 
of a system of two charged surfaces in the presence of 
(polyvalent) counterions and simple monovalent salt ions. 
We proposed two different ensembles for the polyvalent 
counterions, i.e., the canonical ensemble, where the num- 
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FIG. 10: (Color online) Rescaled pressure of the grand- 
canonical system with x = 1 when equilibrated initially at 
the inter-surface half-distance 2o=6 and fixed k — 0.5. For 
slow changes in the distance between the surfaces the pressure 
follows the grand-canonical curve (solid black line). For rapid 
changes in the separation the pressure is expected to follow 
the canonical curve (dashed blue line). 



ber of counterions in the system is fixed and determined 
by the parameter 77, and the grand-canonical ensemble, 
where the number of counterions is regulated by the inter- 
surface separation and via a chemical equilibrium with 
the bulk reservoir where the counterion concentration is 
determined by the parameter x- The grand-canonical en- 
semble does not only quantitatively differ from the canon- 
ical one but also exhibits some very important qualitative 
distinctions that we shall discuss further in this Section. 

As noted before the pressure in the canonical case di- 
verges due to large osmotic contribution from counterions 
in the slit at small inter-surface separations. In contrast, 
the grand-canonical pressure approaches a finite value as 
counterions escape from the slit at small separations. 

Another distinctive feature is that there is a prominent 
hump in the canonical case which makes it possible to 
use a Maxwell construction analysis suggesting a phase 
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coexistence behavior. We find no such behavior in the 
grand-canonical case. 

The canonical description of our model system may be 
relevant for two large thin membranes permeable to sim- 
ple salt (small) ions and impermeable to (large) coun- 
tcrions. It may be assumed that both membranes are 
sealed together at the edges so that the counterions can- 
not escape from the inter-surface region. On the contrary, 
in the grand-canonical ensemble the system is open and 
equilibrates with the bulk reservoir. In practice one must 
of course pay attention to the relaxation time scales since 
polyvalent counterions, being usually bulkier, may typi- 
cally require more time to equilibrate after every change 
of separation than the small salt ions. 

Two equilibration time scales thus unavoidably come 
into play. At every change of separation counterions first 
redistribute in the slit, which happens rather quick. Ad- 
ditionally, they also equilibrate with the counterions in 
the reservoir, but this reservoir equilibration time is usu- 
ally much larger than the first one, especially if the plates 
are very large. For time scales smaller than the reservoir 
equilibration time of counterions to enter the slit, the sys- 
tem would behave as predicted by the canonical variant 
of the theory. 

Figure [TU] shows an example of the pressure-distance 
behavior for a grand-canonical system that was equili- 
brated with the bulk at inter-surface half-distance do = 
6. The amount of counterions in the slit is therefore 
7 ?Gc(5'o)- Thus any change of separation, after the sys- 
tem is left for a sufficiently long time to relax, results 
in the grand-canonical pressure curve as shown by the 
solid line in the figure (Eq. ([52])). On the contrary, if the 
separation changes (starting at the initial value do) are 
fast enough, i.e., for time scales smaller than the counte- 
rion equilibration time with the reservoir, the amount of 
counterions r\ would remain unchanged as 77 = r/cc (Sa- 
lience the system would behave as a canonical system. 
Using the canonical expression (|38p for fixed 77, we can 
calculate the pressure as p(a) — j\ 2 which is shown by 
the dashed line in Fig. [TUJ (note that only the slit re- 
gion behaves canonically and we still need to subtract 
the osmotic bulk contribution of counterions — \x 2 \ this 
term is not present in Eq. (|38[) where no counterions are 
considered in the bulk). When the separation a is in- 
creased in the canonical approach from the equilibrium 
point ao, the osmotic countcrion contribution decreases, 
whereas the outer bulk contribution — \x 2 remains un- 
changed. Therefore, the canonical pressure goes to neg- 
ative values due to the pressure of the counterions from 
the outside and eventually limits to the bulk counterion 
pressure — \x 2 - The deficit of counterions in the slit will 
be compensated by counterions from the bulk only after 
a long enough time, after which the pressure curve would 
coincide with the grand-canonical result. 



B. Finite counterion size 

In order to bring up the principal features of our 
approach, we have assumed here that the counterions 
are point-like. In general, in the SC limit, counterion- 
counterion interactions do not enter the leading-order 
SC theory and thus excluded- volume interactions do not 
affect the leading-order results (reflecting the fact that 
counterions are highly isolated in large correlation holes 
at high couplings) [j, 0, However, counterion size 
effects may become important when additional salt is 
present in the system. 

In practice all ions have finite radii, which, for mono- 
valent salt ion, are typically around 0.1 nm, whereas 
for polyvalent ions, the ion size (radius) R c may be as 
large as 0.3 — 1.0 nm or even larger depending on the 
ion type. Since the (polyvalent) counterions can thus be 
much larger we will give some remarks on the possible 



influence of the countcrion size on our results |4l| . The 
first effect is the hard-core repulsion of counterions from 
the surfaces. This effect is taken into account exactly 
within the leading-order SC theory [1, [39[ via the clos- 
est approach distance between the surfaces which enters 
as an excluded volume in the integrals ([3"U[) and (|53[) . 
This would shift the interaction pressure curves by 2R C 
to the right. On top of that depletion effects due to the 
exclusion of counterions would be visible in the grand- 
canonical case for inter-surface separations below 2R C 
since in that case the counterions cannot enter the inter- 
surface region (42| . 

Note also the electrostatic correction introduced by the 
finite counterion size may be incorporated within the DH 
scheme employed here simply via the well-known charge 
renormalization factor cxp(«;a)/(l + na). For counterion 
radius in the range from 0.2 nm to 0.4 nm and k from 
0.5 nm -1 to 2 nm -1 , this would give na in the range 
0.1 to 1, and thus the renormalized counterion charge 
is varied only by a factor between 1.1 and 1.4. This 
charge valency renormalization therefore leads to only 
small changes in the corresponding value of the coupling 
parameter S (which is proportional to q 3 ) in the range 
1.3 to 3. Another related effect is the failure of the linear 
DH screening in the vicinity of polyvalent counterions 
due to high values of the electrostatic potential [43| . As 
a result, the electrostatic potential in the vicinity of the 
counterions decays faster than predicted by the linear 
screening theory. At larger distances from the counterion 
the electrostatic potential is again DH-likc but with a 
different effective charge value that can be captured by an 
effective valency, q c g. Typically if the countcrion valency 
and radius R c are such that qi^/R c < 10, the effective 
charge is of the order of the bare charge. But above this 
threshold, the effective charge gradually becomes smaller 
than the actual charge and eventually saturates at some 
fixed value q^g [H|,[45|. For typical polyvalent counterion 
size R c ~ 0.6 nm, the effective charge renormalization 
becomes important for a > 9. 

The excluded-volume repulsion between counterions 
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results also from their finite size and leads to structuring 
of the Coulomb fluid at high counterion densities. These 
effects arc however well known and are not distinctive 
for the system at hand. We thus do not venture into the 
details of the hard-core structuring between hard walls 
since it has been amply discussed elsewhere f28j ] . 

There are several other effects that come into play be- 
cause of the ion-ion interactions and correlations such 
as the Bjerrum pair formation which is another manifes- 
tation of the non- linear screening effects (3(| [37| ■ Such 
effects were studied separately from corrections to tradi- 
tional PB and DH approaches. They are expected to be 
important again for qt^/ R c > 10, which falls well outside 
the regime of validity of the dressed counterion theory 
that we consider here. Thus we expect that these finite 
size effects lead to no new physics and no new qualita- 
tive features emerge in the behavior of the system in the 
regime relevant to our study. The details of their possible 
effects will be considered in a separate publication. 



tion becomes larger for smaller screening parameters, k, 
and larger amounts of counterions, rj, in the slit. At 
larger separations the electrostatics is screened out and 
the counterions behave as an ideal gas. Comparison with 
MC simulations reveals that our theory is applicable for 
small enough and large enough separations between the 
surfaces. As expected, the theory also works better at 
larger coupling parameters. 

The grand-canonical case, on the other hand, possesses 
some qualitative differences compared to the canonical 
one. At small separations all the counterions are ejected 
from the inter-surface region and the resulting pressure 
remains finite in contrast to the diverging contribution 
of the counterion osmotic pressure in the canonical case. 
The onset of attraction in the grand-canonical ensemble 
depends on the screening k as well as the bulk concentra- 
tion of counterions \- The MC results confirm that the 
SC theory is relevant for high enough 3 and small surface 
separations as expected from the SC validity criteria. 



IX. CONCLUSIONS 

In this work we have derived a theory to effectively 
solve the model of two charged plane-parallel surfaces 
in the presence of a highly asymmetric electrolyte com- 
posed of monovalent salt and in general polyvalent coun- 
terions. Due to its asymmetric nature we used different 
approaches to describe different components of the elec- 
trolyte solution. The monovalent salt ions were treated 
on a weak-coupling level and the (polyvalent) counteri- 
ons on the strong or weak coupling level. This allowed us 
to formulate a dressed counterion theory that we solved 
for various cases: in the WC and SC limits as well as 
within the canonical and the grand-canonical ensembles 
for the dressed counterions. In the dressed counterion 
theory the salt ions implicitly screen all the interactions 
in the system turning them from a long range into a short 
range screened (DH) form. 

While salt ions were always treated grand-canonically, 
counterions were treated either canonically or grand- 
canonically. The SC dressed counterion theory contains 
only two free parameters, namely, the salt screening pa- 
rameter k and the amount of counterions r\ in the slit in 
the canonical ensemble or \ (corresponding to the bulk 
concentration of counterions) in the grand-canonical en- 
semble. The model itself contains a third free parame- 
ter, namely, the coupling parameter 3 that measures the 
counterion correlation effects and is considered as 3 — > 00 
in the SC dressed counterion theory. The SC dressed 
counterion theory is expected to be applicable at finite 
couplings if the parameter 5 is high enough. In order to 
test the relevancy of the theory we performed also MC 
simulations and compared them with the theory. 

In the canonical case the interaction pressure between 
the surfaces can become attractive at small to interme- 
diate separations due to strong correlations induced by 
polyvalent counterions. The magnitude of the attrac- 
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Appendix A: Simulations 

Simulations of the model system were performed us- 
ing the standard Metropolis Monte-Carlo scheme [46j |. 
appropriately adjusted in the case of a grand-canonical 
system [4?| . when so required for the dressed counteri- 
ons. The salt ions are treated implicitly using the DH 
pair potential (|29p between all other charges. 

In the grand-canonical case, bulk densities were estab- 
lished by separate simulations of a bulk system. In the 
slit geometry, periodic boundary conditions were applied 
in the x — y directions parallel to the charged surfaces. 
The size of the simulation box in these directions, L, 
was always chosen to exceed twice the maximum sepa- 
ration. Furthermore, in all cases we assume L ^> k^ 1 
(in most cases, L was considerably larger than 20k -1 ). 
We nevertheless included a long-range interaction cor- 
rection potential, 3>lr(z), based upon the average coun- 
terion distribution in the slit, (po(z)) [H|]. Specifically, 
the rescaled interaction energy between a counterion at 
a position z and the external (mean-field) distribution is 
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q^LR(z), where 



$ LR (z) = 2^£b^~ 1 / (p (^)>c"^ L2+{2 - 2 ' )2 d^ 



(Al) 

Given the typically large size of the simulation box, such 



a long-range correction is most likely redundant, but 
its implementation on the other hand is computation- 
ally very cheap. The interaction pressure was calculated 
across the mid-plane as well as on the surfaces. The for- 
mer choice leads to slightly better statistical performance 
for hard surfaces. 
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